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Abstract 

We present a numerical method for handling the resolution of a general transport equation for radiative particles, 
aimed at physical problems with a general spherical geometry. Having in mind the computational time difficulties 
encountered in problems such as neutrino transport in astrophysical supernovae, we present a scheme based on full 
spectral methods in 6d spherical coordinates. This approach, known to be suited when the characteristic length of 
the dynamics is much smaller than the domain size, has the potential advantage of a global speedup with respect to 
usual finite difference schemes. An analysis of the properties of the Liouville operator expressed in our coordinates is 
necessary in order to handle correctly the numerical behaviour of the solution. This reflects on a specific (spherical) 
geometry of the computational domain. The numerical tests, performed under several different regimes for the equa- 
tion, prove the robustness of the scheme: their performances also point out to the suitability of such an approach to 
large scale computations involving transport physics for mass less radiative particles. 
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1. Introduction 

Particle transport phenomena are central in modelling systems governed by radiative hydrodynamics, encountered 
very often in astrophysics as well as in plasma physics. A global description of radiative transport involves the 
hyperbolic transport equation, sometimes called Boltzmann equation in the literatur^] This equation describes the 
time evolution of a distribution function F defined on a 6-dimensional phase space. The high dimensionality of this 
equation often prevents its numerical resolution in the most general geometry, due to unaffordable computational 
resources for obtaining physical results in a reasonable CPU time, when using classical techniques. As a result, 
most numerical models for radiative transport physics in several settings either restrict the global geometry of the 
problem (as in [ 18 , 8l fT4l[T9ll . or replace the transport equation by simplified models usually involving the distribution 
moments |fl] [T3] |2T] |6l Q3] . However, in the particular setting of neutrino transport in astrophysical supernovae, the 
fact that a multidimensional transport model for neutrino radiation was needed to reproduce the observed supernova 
explosions has been strongly hinted in recent simulations Ifl4l . As a result, attempts have been made in this direction, 
but they result in very demanding simulations, that are still not able to capture all the needed physics in a general 
geometry |[T7l[T6l . 

Apart from the problem of dimensionality and size of the simulations, transport phenomena very often involve 
physical processes occurring through several orders of magnitude for typical lengths. Once again, this problem arises 
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1 The main difference between the two concepts is that in the transport equation, the collision term only describes interactions between neutri- 
nos/photons and external medium (atoms, nuclei or electrons). (Herein after, we shall refer to neutrinos and photons as "radiative particles", only 
using the term photon or neutrino when it turns out to be necessary). In the Boltzmann case, the collision terms also describe in principle interac- 
tions between the described radiating particles; these are not relevant in our context. Therefore, we will try to avoid the term "Boltzmann equation" 
from now on. 
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in the supernovae neutrinos setting, when comparing the mean free path of a radiative particle in the diffusion regime 
and the typical size of the system iflOl . As long as the two computational problems mentioned above are concerned, it 
is customary to privilege numerical methods with a high order of accuracy [12]. The use of multidimensional spectral 
methods 117113115) seems to be especially adequate. 

Concerning the computational size difficulties, a rule of thumb [7 1 claims in fact that for a given accuracy, spectral 
treatment requires five times less grid points per dimension then the ordinary second order finite difference algorithm. 
A computational factor of 5 6 = 15625 could then be gained in modelling the transport equation. Consequently, 
solving the 6-D transport problem in a reasonable time while using a spectral algorithm and massive parallelisation 
seems possible. In this article, we shall describe the numerical methods that can be used in solving the sole transport 
equation for radiative particles and testing such an economy in computational time. 

From a mathematical point of view and if we follow suitable approximations, the neutrino and photon transport 
equations are analogous in a majority of settings. The abundant results on the photon transport equation will be 
however mainly used here for testing numerically our schem^] 

In this work, the following assumptions will be made: 

1) As mentioned above, we assume the mass of the radiative particles to be zero (which for astrophysical neutrinos 
is a fairly reasonable assumption). 

2) Neutrino and photon polarisation states will not be taken into account, and will be averaged on. 

3) The interacting plasma will be assumed to be at the local thermal equilibrium. This physical oversimplification 
will allow for a much simpler treatment, by enabling to introduce the full thermal equilibrium limit in the equations. 

4) For the simplicity of discussion, possible general relativistic terms will not be taken into account. We believe 
this aspect, though physically important in some astrophysical computations, will in no way change the behaviour of 
the numerical scheme, or the mathematical properties of the studied equations. 

Using (6+1) general spherical coordinates in phase space, we will perform a numerical resolution using spectral 
methods based on Fourier/Chebychev decomposition, depending on the type of coordinates involved. This decompo- 
sition will be performed in the physical space (classical spherical coordinates (r, 9, <p)) as well as in the momentum 
space (energy dependence and angular coordinates for the momentum part). Due to possible discontinuities arising 
on the space variable r, we will also propose an hybrid version of the code in which we use finite differences only 
in this dimension. A tentative conservative full spectral version including also a Chebychev decomposition in r will 
be designed and tested as well. We shall show in different contexts relevant results in 5 dimensions at most, and give 
the corresponding CPU time obtained for every simulation, on single-processor runs of an ordinary computer with a 
clock frequency of 2.5 Ghz. 

The paper is organised as follows: in section 2, we present the general mathematical framework of the transport equa- 
tion, alongside with typical physically motivated source terms for the equation. In section 3 we present two asymptotic 
settings for physical transport of neutrinos and photons, namely the coherent transport case and the Fokker-Planck ap- 
proximation. Those two limiting cases will be used as test problems in the ensuing numerical investigations. Section 
4 will present the derivation of the transport equation in our chosen coordinates, as well as the possible mathemati- 
cal problems that arise with this description; a few solutions will then be proposed to handle resolution in the most 
efficient way. Section 5 presents first numerical tests of the designed code, including time evolution for uniform distri- 
bution and full coherent transport using hybrid discretization. Finally section 6 will present a solution to the spectral 
approach in spherical symmetry, with an emphasis on the issue of conservation of the number of particles. 

2. Transport equation in 6-D 

2.1. Context and definitions 

Let f(x,y,z, fx, p y , p z ,t) be the distribution function in the phase space for a collection of particles, expressed 
in Cartesian-like coordinates. The transport equation will quantify the evolution of this distribution function with 
respect to collision terms, that describe the interaction of the radiative particle with other particle species of a plasma. 



2 The only difference between general transport equation for neutrinos and photons, apart from a difference in cross section expressions, is the 
sign in front of the non linear terms accounting for induced processes. 
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The change in number of particles in the elementary volume of phase space D 3 xD 3 p = dxdy dzdp x dp v dp z is then 
described by 

^/DW/?=CT, (1) 
dt 

where D/dt represents the total derivative and CT includes the collision terms. A reasonable assumption for neu- 
trinos or photons is that between collisions with plasma particles (described in the collision terms), radiating pho- 
tons/neutrinos travel in straight lines with no change in energy. This amounts to the absence of global forces acting 
on the radiating particles. In this context, the non-general relativistic transport equation in Cartesian coordinates takes 
the form: 

I^wfUcT. (2) 
c dt ox' 

where 

to 1 = — , <3-<3 = 1, (3) 

\\p\\ 

c is the light velocity, and CT is a collision term that depends also on /. 

2.2. The collision terms 

Three different types of processes are expressed in collision terms: 

• The rate of spontaneous emission of radiating particles by a particular process on the plasma will be expressed 
as S(v, x); here .x models the spatial position, and v is the radiative particle frequency (or its reduced energy 
E/h, where h is the Planck constant). A general assumption is that the emission is isotropic (i.e. matter itself 
does not have a preferred direction). This of course is only valid if our frame of reference is moving with the 
plasma. 

• Absorption processes are expressed by a cross section cr a (v,x) and are also assumed to be isotropic. 

• Scattering processes at a spatial position x, from a radiating particle scattered from coordinates (to, v) within 
(dto, dv) to (c3', V) within (dto' ,dv') is expressed by the cross section o~d [x, to ■ to , v — > v\ Again, in agreement 
with the previous assumption of isotropy for matter processes, this differential cross section only depends on 
the angle between the incoming and scattered radiating particle momentum, via the simple scalar product to-to' . 

We further assume for simplicity that radiative particles interact with only one species of the plasma, of particle 
mass mo. Using this energy scale we define an auxiliary distribution function 

hv *, t 

F = ( ^ / = 7 2 f, (4) 

where the notation y is coined as the dimensionless energy for the radiative particle. This redefinition of the distribu- 
tion function allows for a slightly simpler notation for the interaction terms, while the left hand side operator of Eq. [2] 
keeps the same form as applied to F. Once all those quantities characterising the interactions with matter are known, 
the collision term is determined by the formula: 



CT = n(r,t){S(v)-o- a (v)F 

+ | dv I du o-diy — > v, at ■ a> )F(v , to )[1 + —, zF(v, to] 
Jo JAn 2v 

r°° , r , c 3 

- I dv \ da> o-div — > v , co ■ u> )F (v, <3)[1 + — -F(v ,to')]}, (5) 
Jo J\k 2v 2 

where n(r, t) is the interacting plasma density. For clarity, no spatial dependence of the chemical composition for the 
plasma is assumed; we finally write here interaction terms that are only proportional to the plasma density n. In the 



right hand side part, and besides the absorption and emission terms, the first integral term models the in-scattered 
neutrinos to coordinates to (c3, v) within (da>, dv). The second integral term models the out-scattered neutrinos, from 
(u>, v) to (to', v'). (Td is the differential scattering kernel for interactions. We have also included in front of the scattering 
integrals the quantum corrections due to induced processes for both types of radiating particle^] Only one term for 
each type of process is represented in an attempt for concision. 



3. Different approximations 

3.1. The coherent scattering 

Consider the very low energy regime for the plasma and the radiating particles; the following assumptions are then 
made: 

1) The plasma is at rest in our frame. 

2) The ratio between the scattered particle energy hv and the scattering target rest mass mnc 2 (be it a lepton or a 
hadron) is very small: (y = hv/moc 2 « l).The velocity of scattering plasma particles will then always be neglected 
in this case. 

Under the above assumptions, we crudely approximate that no energy exchange occurs, meaning that the energy of 
the scattered particle is the same as the impinging one. Therefore the scattering kernels writes: 

o~d(y — * v , & • <3 ) = <Td(0 • & )8(v - V), (6) 

where 5 is the Dirac function. Here we give general expressions for the differential and total cross sections <x rf and cr,, 
for photon and neutrino scattering to electrons and hadrons. The total cross section is defined by: 

cr, - j dv I dto o- d (v — > v , a> ■ a> ). (7) 

Jo J\k 

In the photon/electron case, the coherent scattering approximation leads to the well-known Thomson scattering 
cross sections: 

of = 2r 2 (l+(<S-<5') 2 )c5(v-v), o* h =^r 2 e (8) 

where r e = e 2 /(m e c 2 ) is the classical radius of the electron, m e being its mass. 

In the neutrino/hadron interaction case, the differential cross section is usually reduced to the two leading orders 
in the angular decomposition, in the form: 

o* = -L(A + B(a> ■ a>)) 5(v - v), o* t = A (9) 

47T 

where A and B are constant depending on the nucleons. 

3.2. The Fokker Planck approximation 

For a plasma particle, we denote by a — kT jniQC 2 the thermal energy to mass energy ratio. In this section we 
assume that y « 1 as before, a « 1 and that the plasma, at rest in the laboratory frame, is in local thermodynamic 
equilibrium. In this context, we would like to describe low order energy redistribution in scattering processes. This is 
the setting of the Fokker-Planck approximation^] We obtain then for the photon distribution function F = y 2 f ,1 
Eq.(8.62) : 



dF{yu>) 
dt 



+ oj ■ V ld F{y, (3) = n(x, t)cr Th |-F(y, <3) + J d& [l + (3> ■ & ) 2 ] F(y, 



1287T 2 



I dco I da) [l — £j5 ■ a> + {u> ■ w ) 2 - {a> ■ a> ) 3 J d y \F(y, u> F(y, lo )) > , 

J47T JAK ) 



(10) 



3 The plus sign holds for bosons (photons), the minus sign holds for fermions (neutrinos) 

4 This approximation holds for neutrino nucleon collisions at plasma temperatures < 10" K° and for temperature T < lu'S' and photon energy 
hv < . 1 MeV. It is especially relevant in the context of X-ray astrophysics 
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where we used here a number distribution function, as opposed to the energy distribution function I in [20|. The 
equation is written in the reduced length unit of A c = h/moc , A c = 2.42~ 10 cm being the Compton wavelength. As 
before, some variables dependencies in the distribution function are implicit. Absorption and emission terms are also 
not written here. 



By using the Fokker Planck approximation, one can then replace the integral operator on the energy in the Eq.( 10 1 
by a differential operator, much easier to handle numerically. As mentioned in the introduction, the Fokker-Planck 
limit for transport will be considered as a test case in numerical investigations, alongside with the coherent scattering 
limit. 

From the Fokker-Planck equation we can define two typical times: the "isotropisation time" t, s = l/(ncr r ^c) 
describes the typical evolution of the angular distribution in phase space, whereas the "bosonisation time" tb = 
1 /(ancTThc) will be related to dynamical changes in the energy spectrum of radiative particles. Since a « 1, 
tbo » Ti S holds; consequently, during the evolution, F will undergo an "isotropisation" process in a shorter timescale 
than the energy spectrum of the distribution function F will change significantly. If F is homogeneous and depends 
only on t and y, then Eq.( 10 1 reduces after integration on to and to to the Kompaneet equation IfTTI : 

1 OF 



a 



dt 



ncr T h 



2 dF , 2 

ay — +(y 
ay 



2ay)F+ l -F 7 



0. 



(ID 



If we integrate both sides of the Eq.( 1 1 1 on the dimensionless energy y, we obtain, as expected, an equation which 
expresses the conservation of the number of photons: 

d_ 
dt, 



I 



F(y, t)dy = 0. 



(12) 



The steady state solution of Eq.( 1 1 1 is a Bose distribution: 



F(y) 



2y 2 { 



exp( 



y-fi 



a 



)-l 



(13) 



The factor of 2 in the right hand side being related to photon polarisation, p is an integration constant which physically 
represents the chemical potential. 

The steady state solution is then a Bose distribution and not the usual Planck distribution that describes full thermal 
equilibrium. This is due to the fact that we have omitted the absorption and emission terms, and consequently constrain 
the photon number conservation given by the Eq.(|12), This also justifies the term of "bosonisation" introduced above. 



4. The transport equation in spherical coordinates 

4.1. Definitions and properties 

Starting from the quite general expression for the above equations, we specify now the geometry of our setting, 
as well as the attached chosen system of coordinates. Having in mind transport modelling in astrophysical (stellar) 
settings, the most natural geometry for this type of study is the spherical one. We here choose a set of 6-D spherical 
coordinates related to previously defined phase space vectors (r, to), and described by the variables r, 9, <p, y, 0, <I> as 
in Fig. [1] 

The first three variables are the classical 3D spherical coordinates in physical space; and <I> represent the angular 
dependence in the momentum space, whereas y is a dimensionless measure of the photon (resp. neutrino) energy. 
In this system of coordinates, we can, from the expression in 6-D Cartesian-like coordinates, write the Liouville 
operator to • V = X. using Jacobi matrix products for coordinate changes; one has however to keep in mind that in the 
new coordinate set, the angular variables in the momentum space are defined with respect to physical space angular 
coordinates; this of course slightly complicates the calculation. 

In the end, the operator reads( |20|,and references therein): 

L sp h = cos 



dr 

1 d sin sin O d d cos# d 

sin cos <p - — I sin sin sin <P — - 

86 sin6» dd> 30 sin<9 30 



(14) 




Figure 1: Schematic representation of the 6d spherical coordinate in phase space. Note that a freedom exist for choosing the angular variable <D up 
to a constant. 



so that the general transport equation becomes 

1 dF 

-—+£ sph F = CT. (15) 
c at 

Let us note that it is also possible (and useful) to write this equation in a conservative form: see the Appendix A 
for a derivation of it. 

In order to express integrals in source terms, the expression of the vector to in the new system of coordinates is 
now required. We provide the Cartesian components o) x , 0) y , u> z of the vector to as function of 9, <p, and O: 

o) x = cos sin cos + sin cos <£ cos 9 cos tf> - sin © sin <E> sin tf> (16) 

Cj y = cos sin 9 sin <p + sin cos <I> cos 9 sin <p + sin © sin <I> cos <p (17) 

o) z = cos cos 9 - sin cos O sin 9. (18) 

The following properties hold[^] 

0) 2 x + t^ + a)l = l (19) 

-Csph W. v = 0, t sp h = -Csph 0>z = 0- (2°) 

Therefore, if the distribution function F depends only on Co, we have 

£ sph F(OJ x , 0)y, «;) = 0. (21) 



We shall finish this section by noticing that some terms of the Liouville operator £. sp h given by the Eq.( 14 1 are 
singular for r = and 9 = 0, n. Since the operator is itself regular, these terms correspond to coordinate singularities 
that shall cancel each other in the computation. We shall give an example of such cancellations in our case. Consider 



5 In the Cartesian framework, the identities given by Eq. |20| are quite trivial: consider the Liouville operator £. car t in Cartesian coordinate and 
Cartesian components: 

r - i dF 
■Lean — v 

d.Xj 

For F = p x , F = p y or F = p z the above identities are fulfilled. This obviously holds then for any generic system of coordinates. Numerically, the 
relations in Eq.|20j,or Eq. |21| can be used to assess the numerical accuracy of our resolution. 
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a spherical shell in physical space, for which R\ < r < R2 and R\ > 0. Only singularity issues in 6 — Q,n are then to 
consider. We first write a polynomial decomposition of the distribution in Cartesian-like coordinates: 



F(x,y,z,o x ,co x ,a> z ,t)= J] C(t) i j kAB cx i y J z k aj A ,^co': . (22) 

i,j,k,A,B,C 



In 6-D spherical coordinates, the singular terms in the Liouville operator given by the Eq.( 14 1 are 



sin ^ / d d , 

cos cD — - cos 6— . (23) 



sinf? \d(f> SO 
In the above polynomial decomposition, we encounter two cases: 

- For terms associated with coefficients of type C(0i/#)0 (no dependence on u>), a spherical decomposition in 
(r, 8, <p) will lead to ^-dependent terms being factored by sin{0). 

For terms containing powers of a) x , <o y , o) z , their expression in Eqs.( 16|17|18 1 ensures us that compensation will 



occur when the operator in Eq.(23 1 is applied. 

The spectral representation of the considered fields is able to handle directly the specifics of the decomposition 
(see 12 13 for similar examples). 

4.2. A simplified 2 -dimensional case: The = n/2 discontinuity problem 

We illustrate the prominent difficulties encountered in the analysis of this equation with a problem restricted to 
a spherically symmetric shell (Ri < r < R2) and with only coherent scattering allowed. The solution F for the 
distribution function will then only depend on the three variables t, r, 0. We focus here on analyticity issues and the 
problem of boundary conditions. Under the above hypotheses, the transport equation simplifies to: 



IdF dF sin0 3F 

— -w + cos©— — 

c ot or r o© 

+n e (r)[o- 10t F(t,r,&)- | & dif (cos®cos®') F(t,r,®')sin®' d®'\ = 0, (24) 



L &di 



where n e (r) is a plasma density, cr t0 , and Cdif are respectively the total and differential cross section, and integration on 
the momentum angle €>' has already been performed. In order to perform a very simple analysis, we now artificially 
split the differential operator acting on F, so that we retrieve two advection equations. The radial advection part reads: 

lf + cos©f=0. (25) 
c ot Or 

This is a first order equation, associated to an evolution with velocity V = c cos ©. It propagates from the inner region 
of the shell to the outer one if < < n/2 (cos 9 > 0). On the contrary, it propagates from the outer region to the 
inner one if cos © < 0. Consequently, in our geometrical setting, an inner boundary condition at r — R\ has to be 
imposed for < | (incoming flux) and an outer condition at r — R2 for > n/2 (re-entering flux). 
If we now consider the second advection term 

IJF^ineSF 
c dt r 3® 

the analysis is here simpler: propagation occurs always from © = n to = in the computational domain. However, 
the vanishing of sin at = n shows a degenerate behaviour at this point: no advection in © occurs, therefore no 
boundary treatment is needed. 



Coming back to the full Eq.(24 1, it is now expected that regularity issues in the numerical solution will arisqjacross 
the surface = | , due to different radial advective directions on both sides. For example, a boundary condition value 
for F can be freely set to 

F(t,R u ®<-)=A,AeM, (27) 



We describe a function F as regular if it is of class C p with p large enough to have a fast convergence in the spectral expansion. 
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Figure 2: 2d representation of the computational domain. 



whereas the values F(t,R\,® > |) are advected from the computational domain and therefore uncontrolled. To 
overcome the numerical problems associated with this behaviour, we split our computational domain (here, a spherical 
shell) into two angular domains Dj(r e [R U R 2 ],& e [0, §]) and D 2 (r e [R U R 2 ],& e [f ,n]) (see FigM. To ensure 
particle number conservation across the two domains, we must enforce continuity of the flux on = |. This provides 
us with an incoming boundary condition in © to impose for the solution F in D\, 

5. Numerical tests 

We present below specific tests related to the spectral resolution of the (homogeneous or not) transport equation. 
As outlined above, our computational grid covers a physical shell R\ < r < R 2 split into two domains D\ [0 < < |] 
and D 2 [^ < & < n]. A typical value for our domain size is ^ =5. Unless otherwise stated, spectral decompositions 
are performed using a Chebychev representation on the r,0 and y direction, whereas a Fourier decomposition is 
performed for the remaining angular dependencies 0. 

Time evolution is handled using an explicit second-order Adams-Bashforth method and is performed in the co- 
efficient space, using Tau methods for the imposition of boundaries; however, we also experimented with implicit 
resolution for the coherent scattering problem in spherical symmetry (using Tau-like methods for operator inversion 
in the configuration space) and did not run into any particular issue related to this approach. 

5.1. Time evolution of a uniform distribution 

We assume our domain to be filled by a uniform plasma of constant density n = no, which at first is interacting with 
our radiating particles only through coherent scattering. Absorption and emission are disabled (which ensures particle 
number conservation during the computation) and we start with the artificial initial condition for the distribution 
function: 

1 \ 4 „ _ n 



F(0, <f>, 0, t = 0) = 1 1 + 2cj x + co y + -ljA , < < 



F{0, 4>, 0, t = 0) = 0, 0<©>^. (28) 



Taking advantage of the properties of the Liouville operator described in Eq. 20 we know that at any time of the 
computation, £ sp hF = 0. Monitoring the numerical validity of this property is another way to assess accuracy of our 
approach. 
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Time 

Figure 3: Relative particle number conservation for the coherent time evolution of a uniform distribution 




Time 



Figure 4: L norm of the non-isotropic terms in the spectral decomposition of the distribution function: this is computed as the ratio between the 
coefficient of the constant (isotropic) term in the decomposition, rescaled to the sum of all other coefficients. This serves as a reliable marker for 
the isotropisation process in the evolution. 
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Fig[3]presents particle conservation for this setting over time. The slow drift we encounter only occurs at the level 
of computer roundoff. We consistently obtain a relative error in particle number count smaller than 5.10~ 15 in double 
precision, on timescales much larger than the dynamical timescale of the simulation. The isotropisation process of 
the distribution function due to coherent scattering is also displayed on Fig|4] For those results, the number of points 
used is (N r , Ng, N$,N@,N<t,) = (33, 17, 16, 25, 16). A resolution time step takes about 20 seconds in CPU time. 

Using the same initial spatial profile for the distribution, we now allow for energy dependence and non-coherent 



scattering by implementing the energy-dependent source terms set in Eq. ( 10 1. The initial energy distribution is set to 



be a black body one, at a temperature half the one of the plasma (kT /m e P~- .01 in our units). Conservation of the 
number of photons ensures that F will approach a Bose distribution (see Eq. ( 13 i) over time. Using N y = 33 points in 
the energy dimension, a computational time step takes about 33s. 

In Fig. [5] we can observe the transition made from the initial energy distribution to the final one, and appreciate 
the possible observation of a low-energy condensation that is accessible even with a very limited number of points. It 
is obvious that a specific treatment of the low energy regime (by allocating a specific spectral decomposition domain 
to this region, and increasing the degree of spectral decomposition) would be necessary to study such an effect quan- 
titatively ; however, the goal of this work is only to convince oneself that such study is, indeed, possible with limited 
computational resources. 



t-tj 




Figure 5: Evolution of the energy distribution of F, showing the transition to a Bose distribution and a low-energy condensation. The middle and 
right panel display this energy distribution at four arbitrary consecutive times. No smoothing of the data has been performed. 



5.2. 5d coherent transport in a shell 

We consider a spherical shell enclosing a black body radiating particles (photons/neutrinos) through its outer 
surface R = R\. In the computed shell domain resides a plasma with the following arbitrary density: 

n(r,6,(p) = n [l + Mr sin6» cos (p cos 0)] [1 -(1 - r/R^/il - r/R 2 )] s . (29) 

This plasma triggers coherent scattering, but again emission and absorption processes are disabled for simplicity: 
we only want to monitor the behaviour of a transport process. The shell is initially free of any radiating particles, 
and the central object emits continuously a particle flux following a Lambert law; this leads to the inner boundary 
condition for F: 

F(R\ , 0, <p, ©, cD) = F cos 0. (30) 

This problem will be treated spectrally, except for the radial direction r where we use a simple first order finite- 
difference scheme. The reason behind it is a better treatment of the discontinuity and a reduction of the overshooting 
in this direction that inevitably appears. Grid point numbers are (N r , Ng, N,p, N®, N$>) = (300, 17, 16,25, 16), and a 
time step is around 216s wall clock time, again on a single core. Fig [6] presents distribution function profiles at 
different time steps in the case of an optically thin regime (optical depth with respect to the coherent scattering is set 
to zero), or an optically thick one (in our arbitrary units, the optical depth is set to 5). We are able to represent without 
a problem the beaming effect occurring during the 5d transport, which is also coupled to a large attenuation in the 
second case (the loss of luminosity by 5 orders of magnitudes on a short distance is not altering the code precision, 
as can also be seen by a check on the particle number conservation). It is obvious that in the transparent case, an 
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Figure 6: Profiles of the distribution function F along the directions r and © at time t = 4, 2 in our units. The upper left panel represents angular 
profiles assuming zero optical depth for coherent scattering, whereas the upper right one corresponds to an optical depth of 5 in our units. The 
lower panels are radial profiles at different angles for the same time, and for an optical depth of 5. The boundary values correspond to the thermal 
emission of a Lambertian object (a perfect black body here), and a non-re entrant condition for the outer radius. 



excessive beaming will eventually lead to resolution issues in the angular directions ; this issue can be cured again (at 
least locally) by treating a low © region separately (domain decomposition for the spectral treatment) and assuming a 
better resolution at small angles. It is obvious that an open free streaming region cannot be handled by our approach 
in a clean way : one would have then to resort to less sophisticated descriptions of the particle flux, and match to the 
exact solver. We observe a clear discontinuity of F at the edge r — R\ of the domain, assuming a non-zero optical 
depth; it is a consequence of our rather abrupt assertion for the radiation source to be a pure Lambertian object. A 
more sophisticated approach for the source would allow to get rid of such a feature, although this computation proves 
that the code behaves well even in ill-posed settings. 



6. Enhanced spectral treatment in 2-D case: The conservative formulation 

We present finally a spherically symmetric version of an algorithm for a spectral treatment, amounting to the 
2-D case (F(r, 0, t) for the distribution function, and restricted to coherent Thompson scattering interactions. This 
approach is useful to show how the a prospective full spectral treatment should be handled, and how its inherent diffi- 
culties can be overcome. As opposed to what we did previously, we shall now use the conservative form explicated in 
Eq.( A. 2 1) for the numerical representation of the distribution function: We introduce the new function (see Appendix 
A) 

F(r,®,t) = r 2 sin0F(r,0,f) (31) 
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Using the previously defined value for the total Thompson cross section cr T "' in photon scattering, the 2-D conservative 
form of the transfer equation reads, after integration on the O angl^J 



8F 8F 1 dF . 

— — h cos 0— sin 0- h cos &F = n(r)cr 

ot or r \ 30 



Tot 



-F(r, 0, t) + - sin | 1 1 + (cos cos 0') + -(sin sin 0') ) F(r, 0', t)d® 



(32) 



We shall consider the same boundary problem that the one presented in the 5-D hybrid case, namely F(r, 0, t — 0) = 
as initial value, F(Ri , 0, t) = cos sin for < < n/2 and F(R2, 0, t) for n/2 < < n. The time evolution will 
lead to a discontinuous solution in the radial direction. Spectral methods are not suited to handle this kind of problem. 
In order to show that, consider the simple advection equation 

3<t> <90 

^+C— =0, (33) 
at dr 

where C > is a constant and R[ < r < R2, with the initial data <S){r, 0) = and the boundary condition <$>(R\, t) = 1. 
This problem is clearly analogous to ours, although simpler; it is also well known that the solution is an Heaviside 
function cD(r, t) = ®(r- Cf) propagating in our setting from the inner radius R\ to the outer radius R2 at velocity C. We 
expect to obtain a solution to our problem with similar properties. As we have already said, spectral methods are in 
general not well suited to treat discontinuous solutions, except if some algorithm is used to smear out the solution Q. 
In particular, it is possible to introduce viscosity in the spectral scheme, which can be partly treated in the coefficient 
space |2|. However, such a scheme would have severe effects on conservation laws in our case. In what follows, we 
shall show an algorithm which attempts to overcome such difficulties. 

A classical first order implicit time discretisation to the simple advection problem above (the so-called Euler 
method) leads to 

0)J +1 (r) = (t> J (r) - CAt— — , (34) 
or 

where O-' is the value of the solution at time t = jAt. and Af is elementary time interval. We solve the above equation 
by making an expansion in Chebyshev polynomials, imposing boundary values using a classical Tau approach [7 1 and 
with different values of the parameter Af|j 

Fig.[6]shows the numerical and analytical solutions obtained with N r = 65 points in the propagation direction and 
At - 3(R2-R\)/(C Af 2 ) (this corresponds to 3 times the maximal value satisfying the stability Courant condition for an 
explicit numerical scheme). As expected, we observe strong oscillations in the solution due to the Gibbs phenomenon 
occurring at the solution discontinuity. Note that the propagation velocity for the solution seems however to be 
empirically correct. 

Numerical and analytical solutions using the much bigger value At = 24(^2 - R\)/(CNf) are displayed in Fig. [6] 
While oscillations have disappeared, the numerical solution is spread out (the numerical propagation velocity being 
still correct). We have found experimentally that a value of about 

AV=^2-«i) (35) 

gives the best results, as shown in Fig. [6] The above value seems quite independent on the number of spectral radial 
points: It holds for 17 < N r < 257 in this particular problem. On the contrary, we have observed that a second order 
scheme of the type: 

CAtO& +1 : CAtd& 

+ — — =<!>;-— -3- (36) 

2 or 2 or 



7 This formulation allows to us to check the conservation law term by term: After an angular integration on O, the right-hand side and the second 
term on the left-hand side of the Eq. |32| vanish. The integrated first term on the left expresses then exactly the balance in radial flux. 

8 The matrix of differential operator in the Chebychev basis can be reduced easily thanks a linear combination of the lines, to a a tridiagonal 
matrix. This leads to a considerable speedup of the algorithm. 
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leads to an incorrect propagation velocity. 

In order to compare the results obtained with the first order finite difference scheme and the spectral one, we plot the 
values for the derivative 3 r <t> (see Fig. |6j. We recover similar results for the amplitude of derivatives when the grid 
point number ratio between the finite differences scheme and the spectral one is about five (330 points versus 65). The 
size ratio mentioned in the introduction seems then to be verified. Finally, Fig.|8]shows the error on the conservation 
of particles for this scheme, namely the verification of the identity: 



f 2 <D(r, i 

jRi 



j t | <lHr. tulr + C ( <!»( R\.i) - <!>( A' : ./)) = (). (37) 



The above simple and fast algorithm seems indeed to be able to handle correctly the discontinuities in the solution, 
keeping conservative features at the same time. It is tempting to apply it in solving the Eq.([32|. Let us write this 
equation in the following effective way: 

dF j+l 

F[ +l - cos® k -±-At = P[ + AtS J k +l/2 (38) 
where Ft = F(r,® k ,tj), 0^ is the value of the discretised variable 0, and S k contains all differential terms on F, 



taken on the same values, appearing in Eq. (32i. The coefficient in front of d,F depends on the variable 0. As a 
consequence, we cannot a priori define a consistent optimal value At opr for every value of as in Eq.(35 I. Moreover, 



we want to be free in choosing the value of At, which will in general be constrained by a Courant stability condition, 
either related to the transport equation itself or an adjacent hydrodynamic scheme. We proceed then in the following 



way: consider first an time explicit version of Eq. ( 38 1: 



dF J 

= F{ + cos©i— -At + AtSi. (39) 
K or K 

This can be viewed as a set of N@ equations, where terms are evaluated for each value of the discretized angle 0^. 
We define for each of those angles an optimal time step: 

A4,= JS w2 (*2-*i). (40) 
p cos{® k )Nf 



If At k opt > At, we compute the variation 



G[{t) = F k { tj + At opt ) - F k ( tj ) = A**„ cos 0, . (41) 



A simple linear interpolation is then performed to obtain the updated value for F J k + : 



F(r, ® k , tj + At) = 1 = F{ + G{ cos t ^- + A,S[. (42) 



At 

l opI 

The case At > A t opt is treated by introducing an intermediate time interval 



tl = AtlK<At k opl , (43) 

where K is the smallest integer satisfying the above relation, and performing the numerical integration K times per 
global time step. 



In the implicit setting of Eq. ( 38 1 and with At^ pt > At, we slightly correct the previous scheme by defining: 



G J ke (tj) = At k opt cos A . + e(At k opl )\ (44) 

where the real parameter e is tuned in the algorithm so that the partial update 
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^r=^ + <cos©*-AL (45 ) 

a l opt 

satisfies exactly a radial flux balance for particle number. The update is completed by the implicit first-order step 

P^ = P^* + m{ + \ (46) 



The case A t > A t opt is also performed by splitting time updates as in Eq. (43 1. 

In solving the Eq.([32]l in the domain represented in Fig. [2} strong oscillations due to discontinuities may appear 
near the edges of the interval, if the plasma density there does not vanish (see section 5.2 and Fig.[6]i. A spectral 
resolution is in principle not able to handle these oscillations. To overcome this problem, We have split the interval 
in 3 sub-intervals, two of them being close to the radial edges of the main interval, and each outer sub-interval having 
a width of 0.1 of the main interval. The solution in the sub intervals is computed with a finite difference scheme, 
using 30 grid points. The solution in the largest central interval is computed with the spectral scheme presented just 
above, using 65 spectral points and again performing an expansion on a Chebychev polynomial basis. A Chebychev 
polynomial expansion is also, as before, used to treat the © dependence. 

Results on Fig.|9]show the conservation of the number of particles using this approach, with the settings of Section 
5.2. This shows the validity of our scheme in the bulk and at domain boundaries, and the accuracy of the conservative 
formulation in this 2D example. 

Once known F — r 2 sin &F, the distribution function F is recovered easily by manipulating the coefficients in the 
spectral decomposition; the division by sin is nicely handled in the coefficient space, whereas the division by r is 
performed in the configuration space. 

In conclusion, treating the transport equation with a fully spectral code is not straightforward. In the above (hybrid) 
example in the radial direction, the simulation requires a total of 125 grid points in r. As shown previously, in order 
to obtain the same accuracy with a first order finite scheme, a rough number of 330 grid points would be necessary. 
This leads to a size ratio of 2.64, two times less than previously expected. 

It is possible that the advantage of a spectral scheme reduces with more sophisticated higher order finite differences 
schemes. However, to match the performances of the presented approach, such a scheme should be of order 2 or more, 
exhibit weak diffusivity and show no oscillations due to discontinuities. 

7. Conclusion 

The aim of this numerical work was to assess a "proof of principle" for the treatment of the full transport equation 
in 6D spherical coordinates in a single core processor, in reasonable physical and computational situations, and by 
means of the use of spectral methods in phase space. A particular setting of the computational grid is necessary for 
treating singular behaviour of some terms in the Liouville operator. Meaningful numerical results are obtained in a 
very reasonable computational time, the most time consuming operation being the computation of the Liouville oper- 
ator. For problems where Fokker-Planck-like approximations can not be used, it is possible that the most consuming 
computation would be related to the collision term, in which the thermal distribution of the plasma has to be taken in 



to account. (See Eq.( 14 1. We believe that spectral methods are suited to build an efficient algorithm for the treatment 
this problem. We have also seen that the use of a fully spectral scheme in treating the advection term can turn out to 
be useful if reduction of the CPU. time is a priority. We believe that by using fairly reasonable parallel computation 
on a small-scale cluster, one would be able to perform multiple runs in physically relevant 6-dimensional settings 
and in a really quick fashion. Although we are aware of the fact that several ingredients are still to be added to the 
transport description to use it in a physically relevant radiation hydrodynamics code, our results support the fact that 
no fundamental technical difficulty should arise in tackling those issues. 
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Appendix A. Appendix: Particle conservation in the transport equation 

We will concentrate here on the pure coherent scattering case of transport equation for photons, which writes: 

- ^ + Z sph F + cr,F - f cr d {u> ■ a>)s,m® d<$> = 0. (A.l) 
c dt J 4n 

The above equation can be written, after multiplication by the element volume of the space of phase r 2 sin 6 sin as 

1 d d 

- — (Fr 2 sin 6 sin ©) + — (Fr 2 sin 6 sin cos ©) (A.2) 

c at v ' or v ' 

+ — (Fr sin 2 0sin6» cos U>) (Fr sin 2 sin O] ( A3) 



86 v > d® 



+ — ( F r sin 2 © sin <&) ( Fr sin 2 cos 6 sin <$>) (A.4) 



+ r 2 sin0sin#|o-,F- I cr d (t3 ■ cS)F sin c/0 c/O I = 0. (A.5) 



ie&me(o-,F - 

After an integration on r, 8, <p, 0, O, and provided the detailed balance condition EU1 

0"d( 7 -» 7. <3 -> <3) = cr d { y -> y , <3 -» w ) (A.6) 



holds, we obtain 



where 



N 



^+Jr 2 - Jr, = 0, (A.7) 
f r 2 dr f sm6d6d<f> f sin@d<S>F (A.8) 

JR, J4n J An 



is the number of particles and 



J Rl =R 2 I sm0d6d<p I sin©cos0c/0t/(l)F(f,^2,6i,^,0,a)) (A.9) 

J4tt J4n 

is the flux of ingoing (outgoing) particles into the surface r — R2 of the spherical shell Ri < r < R2- The same 



definition holds for Jr 1 with respect to the radius R\. This is a conservative form of the transport equation (A.l 



In the general case when energy dependence is taken in to account, the photon number conservation is obtained 
after integration on the energy y = hv/mc. If induced scattering processes are also considered, The conservation 
equation Eq (A.7 1 is obtained thanks to the detailed balance condition (A.6 1. 
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Figure 7: Resolution of the advection problem for a shock profile using spectral methods. The left panels show the radial profile for different 
choices of the time step, along with the analytical solution (interpolated to the spectral grid) and a finite-difference solution (lower left panel). In 
the spectral solver, the number of grid points is N r = 65, whereas in the finite-difference solver, the choice N r = 325 is made. The right panel 
shows a comparison for the values of the function derivative in the spectral cases and the finite difference one. The amplitude of the derivative is a 
good indicator of the diffusivity of the numerical scheme. Note that the derivative of spectral solution is larger than the one of the finite difference 
scheme. The value C=l is here chosen. 
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Figure 8: Relative variation of particle number at a given time step, 
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Figure 9: Relative variation of particle number at a given time step, for the hybrid resolution of the 2D transport equation J32J. The number of 
spectral points is N r = 65., with two finite difference zones of 30 radial grid points at the edges. The error is mainly due to the differential term on 
0. It decreases exponentially when the number of grid points in @ increases. 
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